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Summary. We propose a principle of consistency between different hierarchical 
levels of biological systems. Given a consistency between molecule replication and 
cell reproduction, universal statistical laws on cellular chemical abundances are de- 
rived and confirmed experimentally. They include a power law distribution of gene 
PH ' expressions, a lognormal distribution of cellular chemical abundances over cells, 

\^ \ and embedding of the power law into the network connectivity distribution. Second, 

given a consistency between genotype and phenotype, a general relationship between 
phenotype fluctuations by genetic variation and isogenic phenotypic fluctuation by 



r^ , developmental noise is derived. Third, we discuss the chaos mechanism for stem cell 



differentiation with autonomous regulation, resulting froin a consistency between 
cell reproduction and growth of the cell ensemble. 



^^ • 1 Introduction 

o' 

oo ' Biological systems generally form a hierarchy. Ecological systems consist of a 

population of organisms, an organism consists of an ensemble of cells, and a 
t'^^- ' cell consists of interacting biomolecules. Of course, such hierarchical structures 

^D [ also exist in nonliving systems. Then, is there some characteristic property 

underlying biological hierarchical systems? In a hierarchical system, the de- 
scription of units at a lower level and a description of how they come together 
may lead to an understanding of the upper level. However, this bottom-up 
picture may not be sufficient for a complex biological system, since each unit 
^ • at a given lower level is not rigidly determined but can change in adaptation 

5^ I to feedback from a higher level. 

As an example, consider a cell in a multicellular organism, which has in- 
ternal degrees of freedom and can change its chemical composition or gene 
expression patterns. (This is in strong contrast with an electron functioning 
as a hierarchical unit in a physical system). Through interactions with other 
cells, the characteristics of a cell are changed through the process of cell dif- 
ferentiation. A cell in isolation and a cell in a community sometimes exhibit 
different characteristics, since the importance of cell-cell interactions is so 



significant. A cell, and a tissue as an ensemble of cells, mutually determine 
their character. In other words, the character of a unit (e.g., a cell) is deter- 
mined not independently but is changed dynamically by an ensemble of the 
units (see Fig.l). Such dynamic circulation is an essential characteristic of a 
complex biological system [1] [2]: genes encoded in the DNA control macro- 
scopic phenotypes in an organism, while competition between phenotypes at 
the population level determines the expression of genes. 

This interdependence between hierarchy levels has been studied in statis- 
tical physics, in particular, in collective phenomena. Self-consistent solutions 
or approximations are usually adopted in studying cooperative phenomena, 
where stationary, consistent relationships between microscopic elements and 
their mean (collective) field are generated. Although statistical physics is im- 
portant for studying complex biological systems, an essential factor in bio- 
logical systems is not addressed in standard statistical physics: A biological 
unit usually has the potential to reproduce. With reproduction, the number 
of units increases, which may change the relationship between levels, since 
the upper level consists of a population of lower-level units, and this num- 
ber changes in time. This may lead to instability in the consistency between 
elements and the mean (collective) field. Despite changes in population size, 
biological systems generally maintain a degree of consistency between levels, 
even though each unit has many degrees of internal freedoms (e.g., a cell has 
a huge variety of molecules). 

For example, in cell reproduction, the duplication of molecules in a cell is 
correlated such that it keeps some synchrony with the reproduction cycle of 
a cell. In the development of a multicellular organism, reproduction of cells is 
correlated so that the growth of each cell does not interfere with the growth 
of an ensemble of cells. 

Complex System 
^ ^__^ higher level 



^j^-' ^-ie>^J ; 

"^v^ lower level y^ 

flexible units ^ feedback 

Fig. 1. Schematic picture for a hierarchical system. In contrast to a simple system, a 
complex system is regulated by feedback from an upper to a lower level. Consistency 
between the levels should be considered. 



Besides the potential for reproduction and internal reproduction, a biolog- 
ical unit often has the potential to evolve, which requires consistency between 



the time scales of evolution and of development of each unit. Phenotypes 
are generated as a result of the developmental process, which is robust both 
against noise in the developmental process and against some genetic muta- 
tions. Although the time scales of development and evolution are different, 
both are robust against noise or mutation, which suggests consistency be- 
tween development and evolutionfSHHIS]. 

Here we propose that understanding "consistency" among levels and its 
consequence is important for understanding a biological system. In question is 
how such consistency between different levels is sustained and whether there 
are resulting universal laws that apply to all biological systems. 

Here we attempt to answer these questions by considering three examples: 
statistical laws representing consistency between molecule replication and cell 
reproduction; general relationships between genetic variation and phenotypic 
fluctuation resulting from consistency between developmental and evolution- 
ary stability; and the general robust cell differentiation process resulting from 
consistency between cell reproduction and multicellular development. 



2 Consistency between cell reproduction and molecule 
replication 

2.1 Reaction netv^rork for cell reproduction 

A cell consists of several replicating molecular species that help in the synthesis 
of new molecules through catalytic reactions. As a result, a cell grows until it 
divides to produce two cells with similar chemical compositions (see Fig. 2). 



Fig. 2. Basic structure of a reproducing cell with internal catalytic chemical reac- 
tions. Cell-cell interaction is discussed in 54. 



These molecule replications must be synchronized; otherwise, the cellu- 
lar chemical composition cannot be maintained, and the reproduction of cells 
with similar compositions will not continue. At the very least, a membrane 



to separate a cell from the environment must be synthesized, and this pro- 
cess must be synchronized with the replication of other internal chemicals. 
It is unclear how such recursive production and chemical diversity can be 
maintained. Perhaps there is some statistical law for a system to sustain such 
reproduction. 

To investigate the intracellular dynamics of replicating cells, we stud- 
ied several cell models with intracellular catalytic reaction networks that 
transform nutrient chemicals into other chemical species[6] (see Fig. 2 for a 
schematic representation) . Within a cell there are a large number of chemical 
species that mutually catalyze a reaction. There are nutrient chemicals that 
are transported from the outside of a cell and are then transformed to other 
molecules through catalytic reactions. When the number of specific (or total) 
molecule species increases beyond a given threshold, the cell divides. 

We studied a variety of models within this class, by adopting a stochas- 
tic simulation for the reaction. By varying the speed of nutrient intake, we 
found that the cell continues reproduction, keeping an approximately stable 
composition of chemicals, where the growth is optimized[6]. The recursive 
production is satisfied either by tuning the diffusion constant D of nutrient 
across the membrane, or by introducing an active transport process [7]. In 
the former case, the fidelity of reproduction, i.e., the similarity of chemical 
components is maximized at D — *■ D(., where Dc is a critical point beyond 
which the production of a cell does not continue. Around D ~ D^, recursive 
production of a cell with chemical diversity is possible. 

2.2 Universal power law in chemical abundances over species 

We investigated the universal statistical characteristics of reproduction state. 
First, we studied the statistics on the abundance of chemicals for a cell un- 
dergoing reproduction with constant chemical compositions. We measured 
the rank-ordered distributions of chemical species by plotting the number of 
molecules rii as a function of their rank as determined by n^. As first reported 
in [5], the distribution displays a power law with an exponent of -1, often 
called Zipf's law. 

In our model, this power law of gene expression is maintained by a hier- 
archical organization of catalytic reactions. Major chemical species are syn- 
thesized through catalysis by less abundant chemicals. The latter chemicals 
are synthesized by chemicals with much less abundance, and this hierarchy of 
catalytic reactions continues until it reaches the most minor chemical species. 
Indeed, with the aid of mean-field analysis in statistical physics, we can show 
the appearance of a power law distribution with an exponent of -1. 

We have confirmed the universality of Zipf's law by examining a variety 
of conditions in a cell model with reaction networks. They include (i) a dis- 
tributed network connectivity such as the scale-free network distribution, (ii) 
distributed reaction coefficients, (in) higher-order catalytic reactions, (iv) sev- 
eral schemes for transporting nutrients including active transport and passive 



diffusion, and (v) several schemes for cell divisions. Our results are constant 
despite modifications to the model. This power law in regard to abundances 
is observed in any cells that achieve recursive production, i.e., consistency 
between molecule replication and cell reproduction. 

Furthermore, this power law is confirmed by measuring gene expressions 
(i.e. , by measuring the abundances of a variety of mRNAs) . For over a hundred 
cell types examined, we confirmed this power law with an exponent of -1 (see 
also Oil]). The rank-ordered distribution of protein abundances in yeast cells, 
based on protein expression analysis shown infTU^, also suggests the power law 
with an exponent close to -1. 

2.3 Universal lognormal distribution of chemical abundances in 
cells 

We have thus far examined the average abundance of each chemical. Be- 
cause the chemical reaction process is stochastic, the number of each type of 
molecule differs between cells. We therefore studied the distribution of each 
molecule number, sampled among cells, to find that the distribution is fitted 
reasonably well by the lognormal distribution, i.e.: 



p. N 1 / (logn, -logn,)% 

P{n,) exp( ), 1 
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where logrii indicates the average of log rii among cells. In other words, the 
distribution of the chemical abundances is fit by the normal (Gaussian) dis- 
tribution, only after the logarithm of the abundances is taken. 

This lognormal distribution holds for the abundances of all chemicals, 
except for a few chemical species that are supplied externally to a cell as 
nutrients, which obey the standard Gaussian distribution jllj. In other words, 
molecules that are reproduced in a cell obey the lognormal distributions. These 
results beg the question as to why the lognormal distribution law generally 
holds. 

For illustration, consider an autocatalytic process where a molecule (or 
a set of molecules) Xm is replicated with the aid of other molecules. Then, 
the growth of the number rimit) of the molecule species Xm is given by 
dnm{t)/dt = Anm{t), with A describing the rates of the reaction processes 
that synthesize the molecule Xm- Clearly, the synthetic reaction process de- 
pends on the number of the molecules involved in the catalytic process. At 
the same time, however, all chemical reaction processes are inevitably ac- 
companied by fluctuations arising from the stochastic collisions of molecules. 
Consequently, the above rate A has fluctuations ri{t) around its temporal av- 
erage a such that dnm{t) / dt = nm{t){a + il{t)), and hence we obtain: 

d\ognjn{t)/dt = a + r]{t). (2) 

In other words, the logarithm of the chemical abundances shows Brownian 
motion around its mean, as long as r]{t) is approximated by random noise. 



Accordingly, one would expect the logarithm of the chemical abundances (i.e., 
molecule numbers) to obey a normal (Gaussian) distribution. 

In our cell model, however, each reaction process is not autocatalytic, but is 
catalytic with the aid of other molecule species, such that the equation is of the 
type drii/dt = arijUi (i ^ i,j)- Thus, the above discussion cannot be directly 
applied. However, a multiplicative reaction process still exists. For example, 
it is the case that species 3 catalyzes the synthesis of species 1 and 2, as given 
by dni/dt — an2n^, and dn2/dt = anin^. Then, d(ni+n2)/dt = 0(711+^2)713 
follows. If the molecule concentration n^ fluctuates, the above argument on 
multiplicative noise is applied, leading to a lognormal distribution. 

Of course, simple catalysis does not generally exist in our random catalytic 
network. We note a cascade reaction hierarchy, which supports the recursive 
production around the critical state D ^ Dc- A portion of possible reaction 
pathways are used dominantly, which organizes a cascade of catalytic reactions 
so that a chemical in the i-th group is catalyzed by the {i + l)-th, and that 
in the {i + l)-th group is catalyzed by the (i + 2)-th, and so forth. A "modu- 
lar structure" with groups of successive catalytic reactions is self-organized in 
the network. The fluctuations are successively multiplied through this cascade, 
i.e., the noise at the (z -I- 2)-th level multiplicatively influences the {i + l)-th 
level, and {i + l)-th level to i-th level, • • •, and so forth. By taking the log- 
arithm of concentrations (i.e., logrim), these successive multiplications are 
transformed into successive additions of random variables, for which the cen- 
tral limit theorem is applied, leading to the Gaussian distribution of log rim . 

Another simple derivation of lognormal-type distributions in the concen- 
trationof chemicals is provided by noting that the change in concentration 
c of some chemical in a cell is given by dc/dt = (Synthesis)- (Decomposition)- 
(Dilution), where the dilution results from the increase in cell volume V ^ given 
by the term {dV/dt)c. As the growth rate dV/dt fluctuates, there appears a 
multiplicative noise term. By computing a stationary probability distribution 
of c (from the Fokker-Planck equation corresponding to the Langevin equa- 
tion), a log-tailed distribution appears (which approximately agrees with a 
lognormal distribution at a tail). 

The lognormal distribution is observed in a variety of models where the 
cells reproduce efficiently. Following the above argument, the distribution re- 
flects a balance between the replication of each molecule and the growth of 
the cell, leading to the equation drim/dt. This distribution is a result of con- 
sistency between replication of molecules and reproduction of a cell. 

2.4 Embedding the abundance power law into network topology 

Next, we investigated the relationship between the network connectivity 
statistics and the abundance statistics. The distributions in the connectiv- 
ity of reaction networks has been studied extensively ^^ [T5] , while the power 
law in chemical abundances discussed here is independent of the network 
structure, as long as the cell satisfies efficient and recursive growth. 



The stability of reproduction as well as the growth speed may differ based 
on the network structure. We studied the evolution of the network by gener- 
ating slightly modified networks and then selecting those that grew faster. We 
prepared n cells with a randomly connected catalytic network with a given 
initial path number. Then, from each of these mother cells, m mutant cells 
are generated by randomly adding one reaction path to the reaction network 
of the parent cell. Then, for each of the n x m cells, reaction dynamics are 
computed, to obtain the growth speed of each cell. Among the cell population, 
n cells with higher growth speeds are selected. Again, from each of these cells, 
m mutants are generated. This mutation-selection cycle is then repeated. 

In the beginning, the parameters are not set at the critical point, such that 
Zipf 's law on abundances is not observed. By using selection experiments to 
choose cells with higher growth speed, however, Zipf's law on abundances is 
achieved within 10 generations. The network structure is still random. When 
we continue selection processes, however, P{k), the network connectivity dis- 
tribution of reaction path numbers fc, obeys a power law with an exponent 
close to -3 [ig. 

This scale-free-type of connectivity distribution emerges in this evolution 
because attachment of paths to the chemicals with larger abundances is pre- 
ferred as a result of the selection process. Note that the power law distribu- 
tion of chemical abundances has already been established through evolution. 
A connection between a reaction path and a more abundant chemical is more 
effective in increasing the growth speed of a cell. A change in the growth speed 
by addition of an outgoing path from a given chemical species likely increases 
with its abundance x, since the degree of change increases in proportion to 
the flux of the reaction. Then, the probability qout{x) to have such an out- 
going path after selection will increase with x, even though the addition of 
the path itself is random. If such probability linearly increases with x, then 
the abundance power law is transformed to the connectivity power law as 
P{kout) oc kout~ ■ Numerically, we found that the probabilities qout{x) are fit 
by qout{x) ex: a;" with a ~ 1/2. In this case, the connectivity distribution is 
given by fc-("+i)/" ^ fc-3[14]. 

It is interesting to note that the power law in abundances emerges first, 
and later, through evolution, it is embedded into the power law in the network 
connectivity. The abundance power law is in regard to the number of proteins, 
while the network is given by genes that ultimately determine whether a 
particular enzyme species that catalyzes a given reaction is present. When 
genes mutate, the network path is changed accordingly. In this sense, what 
we have observed here can be rephrased as "phenotype (metabolic process 
or gene expression pattern) first, and genes later". The power law in the 
abundance of the former is later "assimilated" by the gene network structure, 
as genetic assimilation as proposed by Waddington[15j. 



3 Consistency between genetic variation and phenotypic 
fluctuation 

3.1 Evolutionary fluctuation response relationship 

The result of §2.3 suggests the existence of large phenotypic fluctuations 
among cells with identical genes. In the model, the network and the parame- 
ters are identical, and in the experiment, isogenic bacteria are used. Still, there 
exist large isogenic phenotypic fluctuations. Here we discuss the relevance of 
such fluctuations to evolution, in relation to genotype-phenotype mapping. 

^ ^ — - " " ^ Genotype a ~ ^^ ~ — ^ ^ 



Fig. 3. Schematic representation of evolutionary processes with genotype and phe- 
notype distributions. 



Phcnotypes are determined from genes through the developmental process. 
However, we note that the developmental process from a single genotype does 
not necessarily produce a single phenotype (see Fig. 3). As mentioned, the 
developmental process is generally noisy, and, hence, the phenotype (such as 
abundances of proteins) from isogenic cells (organisms) fluctuates. 

One might think that such isogenic phenotype fluctuations are not related 
to evolution, since phenotypic change without genetic change is not transferred 
to the next generation. However, the degree of fluctuation is determined by 
the gene, and is heritable. Hence, there may be a relationship between fluc- 
tuation and evolution. We have found evolutionary fluctuation-response rela- 
tionships in bacterial evolution in the laboratory, where the evolution speed 
is proportional to the variance of the isogenic phenotypic fluctuation [16J. This 
proportionality was confirmed in a simulation of the reaction network model 
in the last section [19]. The origin of proportionality between the isogenic phe- 
notypic fluctuation and genetic evolution has been discussed in light of the 
fluctuation-response relationship in statistical physics [17l [18]. 

There is an established relationship between evolutionary speed and phe- 
notypic fluctuation. It is the so-called fundamental theorem of natural se- 
lection by Fisher [21] that states that evolution speed is proportional to the 
variance of phenotype due to genetic variation, which is denoted as V^. It is 



the phenotypic variance as a result of the distribution of genes in a popula- 
tion, defined as the fluctuation of variance of average phenotype over individ- 
uals with different genes. In contrast, the evolutionary fluctuation-response 
relationship, proposed here, concerns phenotypic fluctuation of isogenic indi- 
viduals as denoted by Vip. While Vip is defined as variance over clones, i.e., 
individuals with the same genes, Vg is a result of the distribution of genes. 
Hence, the fluctuation-response relationship and the relationship concerning 
Vg by Fisher's theorem are not identical. 

If Vip and Vg are proportional, the two relationships are consistent. Such 
proportionality, however, is not self-evident, as Vip is related to variation 
against the developmental noise and Vg against the mutation. The relationship 
between the two, if it exists, postulates a constraint on genotype-phenotype 
mapping and may create a quantitative formulation of a relationship between 
development and evolution. 

To determine this possible relationship, we again adopted the cell model 
with catalytic reaction networks and applied the genetic algorithm to evolve 
the network to increase a given fitness. Here, the fitness is given by the number 
Hi^ of a given chemical ig, so that reproducing cells with higher n^^ are selected. 
We evolved cells (with recursive production as mentioned in §2), such that the 
concentration of a given chemical increases. We adopted a genetic algorithm 
with a fitness proportional to the concentration of Ci^ . Here, the mutation rate 
is given by the probability that a path in the network is added or deleted at 
each generation. 

As mentioned, for a given network, there are fluctuations in the abundances 
of each chemical. We took the phenotype variable x = log{ni^), since the 
distribution of rii is approximately lognormal, while theoretical studies |16) 
adopt a variable x whose distribution is close to Gaussian. As a measure of 
the phenotypic fluctuations, we computed a variance of x for a network that 
gives peak abundances at each generation. 

To investigate the distributions of phenotypes due to genetic variation, we 
computed the average phenotypes x over isogenic cells. This average pheno- 
type X differs from mutant to mutant, from which its distribution is obtained. 
The variance of x over all mutants computed from this distribution gives Vg , 
while Vip is just the variance oi x from a given single genotype (i.e., network). 

Vg versus Vip is plotted in Fig. 4, which confirms that Vip oc Vg holds for 
each evolutionary process with a fixed mutation rate. As the mutation rate 
/z increases, the slope of Vg/Vip increases, approaching the fine Vg = Vip. 
When fj, is increased at some value ^J-max i mutant populations exhibiting very 
low values of x increase, the distribution becomes flat, and the peak in the 
distribution shifts downwards. Indeed, around /j, « fJ-maxi Vg is the order of 
Vip. For mutation rates beyond ^max, the phenotype distribution is almost 
flat, as shown in Fig. 5, and the value of x after selection cannot increase from 
generation to generation. Evolution no longer progresses. Thus, the evolution 
speed is optimal around ^ ~ fin 
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Fig. 4. The relationship between Kp, the variance of the phenotype x of the isogenic 
ceUs and Vg, the variance of the phenotype x over 10000 mutants of the selected cell. 
Plotted with identical symbols are the relationships over a given evolution course 
with fixed mutation rates as indicated in the figure (based on jTS]). 
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Fig. 5. Distribution of the phenotype x over 10000 mutants, generated with muta- 
tion rates 0.003, 0.01, 0.02, 0.03, and 0.05. When the mutation rate /i approaches 
\Jimax (in this case, pimax is around 0.03), the distribution is fiattened, and the peak 
position shifts downward. 



3.2 Theoretical discussion 



Is it possible to formulate a phenomenological theory to support the relation- 
ship observed in numerical (and partially in in vivo) experiments presented 
in §3.1? 

Here we consider the distribution both in phenotype x and genotype a. 
Through the evolutionary process, the genotype changes from its dominant 
type a = ao, and then the corresponding average phenotype for each genotype 
changes from x = Xq accordingly. To investigate evolution, both with regards 
to the distribution of phenotype and genotype, we introduced a two-variable 
distribution, P(x,a). 



Vip is the variance of x, which can be written as Vip{a) — J{x ~ 
x{a))'^ P (x , a)dx , where x{a) is the average phenotype of a clonal population 
sharing the genotype a, namely x{a) = J P{x, a)xdx. Vg is defined as the 
variance of the average x{a), over genetically heterogeneous individuals and 
is given by V^ = J{x{a)— < x >)'^p{a)da, where p{a) is the distribution of 
genotype a and <x>as the average of x{a) over all genotypes. 

Assuming the Gaussian distribution, the distribution P{x,a) is written as 
follows: 



P{x, a) =iVexp[- ^^ f°^ + C{a ~ ao){x - Xo)/a -^{a- ao)\ (3) 

With iV as a normalization constant. The Gaussian distribution exp(— ^(a — 
ao)^) represents the distribution of genotypes around a = oq, whose variance 
is (in a suitable unit) the mutation rate /Lt. The above equation can then be 
rewritten as: 



„/ N C^ r (X — Xq — C (a — ao))'^ , C^ In/ ■,0^ /,N 

^(^' ") - ^^^p[-^ — 2Mk — - + (2^ - Y,^^- - «°) ]■ (4) 

Our second assumption evolutionary stability states that at each stage 
of the evolutionary process, the distribution has a single peak in (a;, a) space. 
In order for this distribution to have a single peak (i.e., not to be flattened 
along the direction of a), the following condition (besides a > 0) should be 
satisfied: ^ - ^ < 0, i.e., ^ < ^ = ^imax- 

This means that the mutation rate has an upper bound ^max beyond which 
the distribution does not have a peak in the genotype-phenotype space. Be- 
yond this mutation rate, the distribution is extended to very low values of 
X (fitness). This breakdown of the high- fitness phenotype is a kind of error 
catastrophe, if we follow the term by Eigen'20]. This is consistent with the 
observation in Fig. 5. 

We investigated the phenotypic variance due to the genotype distribu- 
tion. First, we considered the average Xa over the distribution P{x, a) for a 
given fixed a, and then considered the distribution of Xa according to the 
distribution p(a), noting that Xa = J xP{x,a)dx — Xq + C{a — ag). For the 
population having identical phenotype x — Xq, the genetic variance in this 
distribution is given by < {5aY >— M. Hence, the phenotype variance from 
this population (of given phenotype x) is written as: 

V^g =< {Xa - Xaof >= C^ fi. (5) 

We note that C / 0, so that the (average) phenotype changes with the change 
of genotype. Noting that a is the phenotypic variance < Sx^ > of isogenic 
individuals, Vip, the inequality /i < a/C'^ is rewritten as: 



V^g < V,p. (6) 

Since the genetic variance in the population < (Sa)'^ > is proportional to 
the mutation rate, the above inequality is simply the threshold for the error 
catastrophe mentioned above. In other words, at the threshold mutation rate 
as Umax, Vig = Vip holds. Then, recalling V^g ex /i, we get Vig = j^Vtp- 

In general, Vg ^ Vig, as Vg is the variance over all current populations. In 
this case, the variance < {5a)^ > is computed over the entire population, and 
is given by ^/(l - ^lCya). Thus, Vg = Vg/{1 ~ Vg/Vp). (Note Vg/Vp < 1). 
If the mutation rate is small, however, Vg w Vig, so that the proportionality 
between Vg and Vip is explained from the above argument. 

To sum up the present formulation of distribution, we have obtained 

(i) Vip > Vig, and (ii) proportionality between Vg and Vip for small mu- 
tation rate cases through a given course of evolution are derived. 

Although these relationships are supported by the above simulation as well 
as by recent studies about gene networks [22], the above formulation is not a 
'derivation'. First, we assumed the existence of two-variable distributions in 
genotype and phenotype P{x, a). As genetic change is not given simply by the 
change of a continuous parameter, it is not a trivial assumption. (For example, 
the genetic change in the cell network model is an addition or a deletion of 
a reaction path, and if expressed as a continuous variable is not self-evident). 
Second, the stability assumption assuring a single peak is expected to be valid 
for gradual evolution. Third, to adopt eq.(3), the existence of error catastrophe 
(to produce mutants with very low fitness values at a large mutation rate) is 
implicitly assumed. With these three assumptions, error catastrophe is shown 
to occur at Vig w Vip. 

In numerical evolution, the error catastrophe is estimated around Vg ~ Vip. 
Here, the phenotype at each generation is within a small range, and the devia- 
tion oiVig from Vg is not so large. Indeed, the estimate of the critical mutation 
rate for the error catastrophe is not accurate enough to distinguish between 
the two. Thus, the above theoretical estimate for the error catastrophe is 
consistent with the numerical result. lJ 

The above results ask why the error catastrophe produces low-fitness phe- 
notypes. Note that the growth of a cell in our model (and in nature) requires 
maintenance of a variety of chemicals through reproduction. Mutants may fail 
to synthesize some chemicals concurrently. When the isogenic phenotypic fluc- 
tuation is large, there is room to search for networks that are robust against 
such mutational change, while for a large mutation rate, the network with the 
highest fitness is not maintained over generations. 

The relationships (i)-(ii) as well as the estimate of error catastrophe are 
also confirmed in a gene network model[22]. We expect these relationships (as 



As Vg 7^ Vig, the inequality between Vip and Vig does not set the bound between 
Vg and Vip. Hence, the bound for heritability as anticipated at the discussion part 
ofll9l is not derived. 



well as the existence of the error catastrophe) to be generally valid for systems 
satisfying the following conditions. 

(A) Fitness is determined through developmental dynamics. 

(B) Developmental dynamics is complex such that its orbit, when deviated 
by noise, may fail to reach the state with the highest fitness. 

(C) There is effective equivalence between mutation and noise in the de- 
velopmental dynamics with regards to phenotype change. 

Condition (A) is straightforward in our model, and condition (B) is satis- 
fied because of the complex reaction dynamics to sustain the growth of a cell. 
Postulate (C) is satisfied since either the phenotypic fluctuation to increase the 
number of catalysts for a synthesis reaction or mutation to add another reac- 
tion term for its synthesis can contribute to the increase in the concentration 
of a given chemical in the same manner. Condition (A)-(B) supports the exis- 
tence of 'error catastrophe'. In several models of evolution, fitness is directly 
represented as a function of genotype. Then, a slight change in genotype does 
not lead to major differences in phenotype. On the other hand, in a system 
with (A)-(B), slight differences in genotype may lead to huge differences in 
phenotype, as the phenotype is determined after temporal integration of the 
developmental dynamics where slight differences in genotype (rule governing 
the dynamics) are accumulated and amplified. 

Note that relationship (ii) shows correlation between the phenotypic 
change by gene and spontaneous phenotypic fluctuation. In other words, the 
degree to which genes can alter phenotypes is presumed in the spontaneous 
phenotypic fluctuation. This suggests consistency between genetic and phe- 
notypic levels (see Fig. 3), as was first discussed by Waddington, as genetic 
assimilation [15] . 



4 Consistency between cell replication and reproduction 
of multicellular organisms 

We briefly discuss cell differentiation, i.e., diversification into a discrete set of 
cell types through development and robustness in the population distribution 
of each cell type through development. 

In cell differentiation in multicellular organisms, there usually exists ir- 
reversibility. Embryonic stem cells have the potential to produce all types of 
cells. Stem cells produced from them have the potential to produce only a lim- 
ited class of cells (for example, a blood stem cell can only generate all types 
of blood cells). The stem cells proliferate or differentiate into other types with 
some probability. Successive differentiation from stem cells (S) terminates at 
determined cells, which can only replicate (see Fig. 6). This cell society must 
then maintain some ratio of population sizes of cell types. Stem cells must 
then regulate the ratio of proliferation and differentiation, depending on the 
population distribution of each cell type. In the example displayed in Fig. 6, 
the ratio of differentiation S —^ A has to be increased when the number of 




Fig. 6. Schematic representation of hierarchical differentiation from a stem cell. 



type A cells is decreased, while the ratio of proliferation S ^ S is increased 
when the number S is decreased. 

We have studied a reaction network model as given in Fig. 2, choosing the 
model so that the chemical concentration changes chaotically [23j. We consid- 
ered cell-cell interactions. As the cell number increases by cell divisions, the 
cells share common resources. Through the concentration change of resources 
in the medium, cells interact with each other. This interaction introduces a 
change in intracellular dynamics. 

Here, cell state is given by the composition of chemicals (i.e., by location in 
the phase space of chemical concentrations) . Each cell type is determined as an 
'attractor' in the intracellular reaction dynamics during cell-cell interaction. 
Initially, the cellular state is a chaotic attractor that gives type 'S'. With cell 
divisions, we found that some cells differentiate from the original cell state 'S' 
to a different type with a distinct chemical composition ('A'), while another 
type ('B') is differentiated later. With further increases in the cell number, 
differentiation from A to Al, A2, and A3 progresses. In summary, cell differ- 
entiation progresses as displayed in Fig. 6. Stem cells differentiate into other 
types stochastically due to the instability caused by cell-cell interactions. 

The differentiation is 'stochastic', arising from chaotic intracellular chemi- 
cal dynamics. The choice of a stem cell either to proliferate or to differentiate 
appears stochastic as far as the cell type is concerned. However, this is not 
due to external fluctuations but is a result of the intracellular state. As this 
state is influenced by cell-cell interaction, the probability of differentiation 
can be regulated according to the population of cell types. As the number 
(fraction) of stem cells increases, the basin for such chaotic attractors touches 
with the attractor itself so that the cells switch to a different cellular state. 
On the other hand, as the fraction of stem cells decreases, the original chaotic 
attractor is stable so that the stem cells can proliferate. With this spontaneous 
tuning in the stability of the stem cell state, the fraction of each cell type is 
regulated. On the other hand, the state of stem cells is 'marginally stable', in 



the sense that its attractor is on the verge of touching with its basin, through 
this regulation. 

Reaction networks ahowing for chaotic dynamics and differentiation are 
not common. We checked a variety of networks to examine the growth speed 
of each cell (i.e., the inverse of the division time) and the ensemble of cells |24). 
First, cells having networks allowing for chaotic dynamics and differentiation 
do not grow quickly at the cellular level. Cells having other networks without 
any oscillatory dynamics in chemical concentrations often divide faster. How- 
ever, as the number of cells grows, the speed of division for such non-oscillating 
cells is drastically reduced, while for cells with chaotic dynamics and differ- 
entiation, the speed is not so reduced. This is because cells of the latter type 
do not strongly compete for resources. In other words, consistency between 
a single cell division and the growth of an ensemble of cells is achieved for 
cells having networks allowing for chaotic dynamics. Hence, the consistency 
between reproduction of each cell type and population growth as an ensemble 
of cells works as a pressure to select networks exhibiting chaotic dynamics. 
Regulation to maintain the proportion of each cell type and the 'marginal sta- 
bility' of the stem cell state are a consequence of the cell system that satisfies 
the consistency. 



5 Conclusion 

Here we reviewed three problems in biology from the viewpoint of 'consistency 
between different levels'. First, as a result of consistency between molecule 
replication and cell reproduction, chemical reaction dynamics are shown to 
be at a critical state, and a power law distribution of chemical abundances 
(gene expression) is derived. The dynamics of the molecule number is shown to 
include a multiplicative noise term, which leads to a lognormal distribution of 
chemical abundances over cells. Through evolution, the power law distribution 
of abundances is embedded into the network topology, leading to a scale-free 
network, demonstrating consistency between reaction dynamics and network 
structure. 

Second, the genotype-phenotype relationship is discussed. We found a 
general relationship between phenotype fluctuations by genetic variation, Vg , 
and isogenic phenotypic fluctuation by developmental noise, Vip, as a result 
of a consistency between genotype and phenotype, because of feedback from 
phenotype to gene through selection. 

As the third topic, we touched upon the chaos mechanism for stem cell dif- 
ferentiation with autonomous regulation, as a result of consistency between 
cell reproduction and growth of the cell ensemble. Although not described 
here, recent studies show the existence of spontaneous adaptation as a re- 
sult of consistency between cellular growth and (stochastic) gene expression 
patterns [25l [26] . We propose that the consistency principle is generally rele- 



vant to understanding reproduction, adaptation, evolution, and development 
in biological systems. 

The authors would like to S. Sawai, M. Tachikawa, K. Fujimoto, and T. 
Yomo for stimulating discussions. 
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